/***
Title
-----

__crsconvert_core__ - Core functionality for coordinate reference system conversion

Syntax
------

> __crsconvert_core__ _xvar_ _yvar_, gen(_prefix_) from(_source_crs_) to(_target_crs_)

Description
-----------

__crsconvert_core__ converts coordinates from one coordinate reference system (CRS) to another.
This command provides the core functionality for coordinate transformation by utilizing
the GeoTools Java libraries.

Options
-------

* __gen(string)__ specifies the prefix for the new variables to be created containing
    the converted coordinates.

* __from(string)__ specifies the source coordinate reference system. This can be:
    - An EPSG code (e.g., "EPSG:4326")
    - A WKT string representing a coordinate reference system
    - A path to a GeoTIFF file (.tif)
    - A path to a Shapefile (.shp)

* __to(string)__ specifies the target coordinate reference system. This accepts the same
    formats as the __from__ option.

Details
-------

The command processes the input coordinates (_xvar_ and _yvar_) and creates two new 
variables with the prefix specified in __gen()__ followed by the original variable names.

For file paths in __from()__ and __to()__, both relative and absolute paths are supported.
If a relative path is provided, it's resolved relative to the current working directory.

The actual coordinate transformation is performed using the GeoTools Java libraries.
When using shapefiles as CRS sources, the files must include their auxiliary components 
(.shx, .dbf, and preferably .prj).

Example
-------

Convert coordinates from WGS84 (EPSG:4326) to UTM Zone 10N (EPSG:32610):
***/

cap program drop crsconvert_core
program define crsconvert_core
version 18

syntax varlist(min=2 max=2 numeric), gen(string) from(string) to(string)

local x: word 1 of `varlist'
local y: word 2 of `varlist'

confirm new var `gen'`x'
confirm new var `gen'`y'

qui gen double `gen'`x' = .
qui gen double `gen'`y' = .

// 处理 from 和 to 参数的路径
local from `from'
local to `to'

// 检查 from 是否是文件路径
if strpos(lower("`from'"), ".tif") | strpos(lower("`from'"), ".shp") {
    removequotes, file(`from')
    local from `r(file)'
    local from = subinstr("`from'", "\", "/", .)
    if !strmatch("`from'", "*:\\*") & !strmatch("`from'", "/*") {
        local from = "`c(pwd)'/`from'"
    }
    local from = subinstr("`from'", "\", "/", .)
}

// 检查 to 是否是文件路径
if strpos(lower("`to'"), ".tif") | strpos(lower("`to'"), ".shp") {
    removequotes, file(`to')
    local to `r(file)'
    local to = subinstr("`to'", "\", "/", .)
    if !strmatch("`to'", "*:\\*") & !strmatch("`to'", "/*") {
        local to = "`c(pwd)'/`to'"
    }
    local to = subinstr("`to'", "\", "/", .)
}

// 调用 Java 方法
java: crsconvert("`x'", "`y'", "`gen'`x'", "`gen'`y'", "`from'", "`to'")

end

cap program drop removequotes
program define removequotes,rclass
version 16
syntax, file(string) 
return local file `file'
end


java:


/cp gt-main-32.0.jar
/cp gt-referencing-32.0.jar
/cp gt-epsg-hsql-32.0.jar
/cp gt-epsg-extension-32.0.jar
/cp gt-geotiff-32.0.jar
/cp gt-coverage-32.0.jar
/cp gt-shapefile-32.0.jar
/cp gt-api-32.0.jar
/cp gt-metadata-32.0.jar

/cp json-simple-1.1.1.jar
/cp commons-lang3-3.15.0.jar
/cp commons-io-2.16.1.jar
/cp jts-core-1.20.0.jar



import org.geotools.api.referencing.crs.CoordinateReferenceSystem;
import org.geotools.api.referencing.operation.MathTransform;
import org.geotools.api.referencing.operation.TransformException;
import org.geotools.geometry.jts.JTS;
import org.geotools.geometry.jts.JTSFactoryFinder;
import org.geotools.referencing.CRS;
import org.locationtech.jts.geom.Coordinate;
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.GeometryFactory;
import org.locationtech.jts.geom.Point;
import com.stata.sfi.Data;
import com.stata.sfi.SFIToolkit;
import org.geotools.gce.geotiff.GeoTiffReader;
import org.geotools.coverage.grid.GridCoverage2D;
import org.geotools.data.shapefile.ShapefileDataStoreFactory;
import org.geotools.data.shapefile.ShapefileDataStore;
import java.util.logging.ConsoleHandler;
import java.util.logging.Handler;
import java.util.logging.Level;
import java.util.logging.Logger;

// Initialization method to replace the static block
void initializeGeoTools() {
    // Disable the JSON-related service loading at startup
    System.setProperty("org.geotools.referencing.forceXY", "true");
    System.setProperty("org.geotools.factory.hideLegacyServiceImplementations", "true");

    // Suppress specific service loader errors
    Logger logger = Logger.getLogger("org.geotools.util.factory");
    logger.setLevel(Level.SEVERE);

    // Suppress INFO level messages from GeoTools
    Logger geoToolsLogger = Logger.getLogger("org.geotools");
    geoToolsLogger.setLevel(Level.WARNING);
    for (Handler handler : geoToolsLogger.getHandlers()) {
        if (handler instanceof ConsoleHandler) {
            handler.setLevel(Level.WARNING);
        }
    }
}


void crsconvert(String x, String y, String newx, String newy, 
                            String sourceCRS, String targetCRS) {

    
    try {
        // Call the initialization method
        initializeGeoTools();

        // Reset CRS cache
        CRS.reset("all");

        // 判断 sourceCRS 和 targetCRS 的类型
        CoordinateReferenceSystem source = parseCRS(sourceCRS);
        CoordinateReferenceSystem target = parseCRS(targetCRS);

        // 打印转换信息
        System.out.println("Converting coordinates from CRS:");
        System.out.println("Source CRS: " + source.toString());
        System.out.println("Target CRS: " + target.toString());

        // 创建转换器
        MathTransform transform = CRS.findMathTransform(source, target, true);
        GeometryFactory geometryFactory = JTSFactoryFinder.getGeometryFactory();

        // 复用坐标对象
        Coordinate coord = new Coordinate();
        Point point = geometryFactory.createPoint(coord);

        // 获取变量索引
        int xIndex = Data.getVarIndex(x);
        int yIndex = Data.getVarIndex(y);
        int newxIndex = Data.getVarIndex(newx);
        int newyIndex = Data.getVarIndex(newy);
        Long TotalObs = Data.getObsTotal();

        // 转换循环
        for (int i = 1; i <= TotalObs; i++) {
            try {
                coord.x = Data.getNum(xIndex, i);
                coord.y = Data.getNum(yIndex, i);

                Geometry transformed = JTS.transform(point, transform);
                Coordinate result = transformed.getCoordinate();

                Data.storeNumFast(newxIndex, i, result.x);
                Data.storeNumFast(newyIndex, i, result.y);
            } catch (Exception e) {
                SFIToolkit.error("Error at obs " + i);
            }
        }
        Data.updateModified();
    } catch (Exception e) {
        SFIToolkit.error("Conversion failed: " + e.getMessage());
    }
}

// 修改后的 parseCRS 方法
CoordinateReferenceSystem parseCRS(String crsInput) throws Exception {
    try {
        // 如果输入是 GeoTIFF 文件
        if (crsInput.toLowerCase().endsWith(".tif")) {
            return readCRSFromGeoTIFF(crsInput);
        } 
        // 如果输入是 Shapefile 文件
        else if (crsInput.toLowerCase().endsWith(".shp")) {
            return readCRSFromShapefile(crsInput);
        } 
        // 如果输入是 EPSG 编码
        else if (crsInput.startsWith("EPSG:")) {
            return CRS.decode(crsInput, true); // 强制使用 XY 顺序
        } 
        // 假设输入是 WKT 字符串
        else {
            return CRS.parseWKT(crsInput);
        }
    } catch (Exception e) {
        throw new Exception("Failed to parse CRS from input: " + crsInput + ". Error: " + e.getMessage(), e);
    }
}


// 从 GeoTIFF 文件读取 CRS
CoordinateReferenceSystem readCRSFromGeoTIFF(String filePath) throws Exception {
    GeoTiffReader reader = null;
    try {
        reader = new GeoTiffReader(new File(filePath));
        return reader.getCoordinateReferenceSystem();
    } catch (Exception e) {
        throw new Exception("Failed to read CRS from GeoTIFF file: " + filePath + ". Error: " + e.getMessage(), e);
    } finally {
        if (reader != null) {
            reader.dispose(); // 确保释放资源
        }
    }
}

// 从 Shapefile 文件读取 CRS
CoordinateReferenceSystem readCRSFromShapefile(String filePath) throws Exception {
    ShapefileDataStore shapefileDataStore = null;
    try {
        File shpFile = new File(filePath);
        if (!shpFile.exists()) {
            throw new Exception("Shapefile does not exist: " + filePath);
        }

        // 检查 Shapefile 的必要组件
        String basePath = filePath.substring(0, filePath.lastIndexOf("."));
        File shxFile = new File(basePath + ".shx");
        File dbfFile = new File(basePath + ".dbf");
        File prjFile = new File(basePath + ".prj");

        /* if (!shxFile.exists() || !dbfFile.exists() || !prjFile.exists()) {
            throw new Exception("Missing required shapefile components (.shx, .dbf, or .prj): " + filePath);
        } */

                    if (!shxFile.exists() || !dbfFile.exists()) {
                System.out.println("Warning: Missing required shapefile components:");
                if (!shxFile.exists()) System.out.println(" - Missing .shx index file");
                if (!dbfFile.exists()) System.out.println(" - Missing .dbf attribute file");
                if (!prjFile.exists()) System.out.println(" - Missing .prj attribute file");
                System.out.println("A complete shapefile requires .shp, .shx, .dbf and .prj files.");
                throw new Exception("Incomplete shapefile: " + filePath);
            }

            ShapefileDataStoreFactory dataStoreFactory = new ShapefileDataStoreFactory();
            Map<String, Object> shpParams = new HashMap<>();
            shpParams.put("url", shpFile.toURI().toURL());
            shapefileDataStore = (ShapefileDataStore) dataStoreFactory.createDataStore(shpParams);

            // Set UTF-8 encoding explicitly
            /* shapefileDataStore.setCharset(java.nio.charset.Charset.forName("UTF-8")); */
        CoordinateReferenceSystem crs = shapefileDataStore.getSchema().getCoordinateReferenceSystem();
        /* System.out.println("PRJ file path: " + prjFile.getAbsolutePath()); */

        /* // 使用 ShapefileDataStore 解析 CRS
        Map<String, Object> params = new HashMap<>();
        params.put("url", shpFile.toURI().toURL());
        ShapefileDataStoreFactory factory = new ShapefileDataStoreFactory();
        shapefileDataStore = (ShapefileDataStore) factory.createDataStore(params);

        if (shapefileDataStore == null) {
            throw new Exception("Failed to create ShapefileDataStore for: " + filePath);
        }

        // 获取 CRS
        CoordinateReferenceSystem crs = shapefileDataStore.getSchema().getCoordinateReferenceSystem();
if (crs == null) {
            throw new Exception("CRS is null for Shapefile: " + filePath);
        }*/

        return crs; 
    } catch (Exception e) {
        throw new Exception("Failed to read CRS from Shapefile: " + filePath + ". Error: " + e.getMessage(), e);
    } finally {
        if (shapefileDataStore != null) {
            shapefileDataStore.dispose(); // 确保释放资源
        }
    }
}

end